
/**
 one can try
 ./md -np 100 -box 5,5,5 -dt 0.0005 -tmax 10 -step verlet -T 10 -therm 3
 ./md -np 108 -box 4.4 -dt 0.0005 -tmax 20 -step verlet -T 150 -therm 3      gaz
 ./md -np 256 -box 2 -dt 0.0005 -tmax 20 -step verlet -T 50 -therm 3		 solid
 **/

/** TODO: - put into .tex calcul cost (N(N-1))...
          - change t_run - t_old to integer counting.
          - algorithm to evaluate forces through cells...
 
	UNITS ADIM: -solid: density=1.2,  temperature=0.5
				-gas:   density=0.05, temperature=1.25
 **/			

#include "./mds_density.h"
#include "../mds/run.h"
#include "../mds/io/xyz.h"
#include "../mds/io/mds.h"
#include "../mds/constants.h"
#include <math.h>

/*******************************************************************************
 *
 *                  Options Callbacks
 *
 ******************************************************************************/
static
int on_option_np( const char *value, void *args )
{
	size_t *num_particles = args;
	*num_particles = strtol(value, NULL, 10);
	if( *num_particles <= 0 )
	{
		fprintf(stderr, "** invalid num_particles\n");
		return -1;
	}
	return 0;
}

static
int on_option_box( const char *value, void *args )
{
	mds_vector *box = args;
	if( mds_box_parse( box, value )  < 0 )
	{
		mds_perror( stderr, "mds_box_parse('%s')", value );
		return -1;
	}
	return 0;
}

static int on_option_dt( const char *value, void *args )
{
	Real *dt = args;
	*dt = (Real)strtod( value, NULL);
	if( *dt <= 0 )
	{
		fprintf( stderr, "** invalid dt=%g\n", *dt );
		return -1;
	}
	return 0;
}

static int on_option_tmax( const char *value, void *args )
{
	Real *tmax = args;
	*tmax = (Real)strtod( value, NULL);
	if( *tmax <= 0 )
	{
		fprintf( stderr, "** invalid tmax=%g\n", *tmax );
		return -1;
	}
	return 0;
}

static int on_option_mass( const char *value, void *args )
{
	mds_mass_type *mass_type = args;

	if( !strcmp( value, "1" ) )
	{
		*mass_type = MDS_MASS_1;
		return 0;
	}
	
	if( !strcmp( value, "alea" ) )
	{
		*mass_type = MDS_MASS_ALEA;
		return 0;
	}
	fprintf( stderr, "** unknown mass '%s'\n", value );
	return -1;
}

static int on_option_units( const char *value, void *args )
{
	mds_units_type *units_type = args;

	if( !strcmp( value, "real" ) )
	{
		*units_type = MDS_UNITS_REAL;
		return 0;
	}
	
	if( !strcmp( value, "adim" ) )
	{
		*units_type = MDS_UNITS_ADIM;
		return 0;
	}
	fprintf( stderr, "** unknown units type '%s'\n", value );
	return -1;
}

static int on_option_step( const char *value, void *args )
{
	mds_step_type *step_type = args;
	if( !strcmp( value, "verlet" ) )
	{
		*step_type = MDS_STEP_VERLET;
		return 0;
	}
	
	if( !strcmp( value, "beeman" ) )
	{
		*step_type = MDS_STEP_BEEMAN;
		return 0;
	}
	
	fprintf( stderr, "** unknown step '%s'\n", value );
	return -1;
}

static int on_option_potential( const char *value, void *args )
{
	mds_potential_type *potential_type = args;
	if( !strcmp( value, "lj" ) )
	{
		*potential_type = MDS_POTENTIAL_LJ;
		return 0;
	}
	
	if( !strcmp( value, "lj-cutoff" ) )
	{
		*potential_type = MDS_POTENTIAL_LJCUTOFF;
		return 0;
	}
	
	fprintf( stderr, "** unknown potential '%s'\n", value );
	return -1;
}

static int on_option_T( const char *value, void *args )
{
	Real *T = args;
	*T = (Real)strtod( value, NULL);
	if( *T < 0 )
	{
		fprintf(stderr, "** invalid T=%g\n", *T);
		return -1;
	}
	return 0;
}

static int on_option_therm( const char *value, void *args )
{
	Real *therm = args;
	*therm = (Real)strtod( value, NULL);
	return 0;
}

static int on_option_f( const char *filename, void *args )
{
	mds_options *options = args;
	fprintf( stderr, "-- appending options from file '%s'\n", filename );
	if( mds_options_load( options, filename, "f") < 0 )
	{
		mds_perror( stderr, "mds_options_load(%s)", filename );
		return -1;
	}
	return 0;
}

static int on_option_noise( const char *value, void *args )
{
	Real *noise = args;
	fprintf( stderr, "-- reading noise='%s'\n", value );
	*noise = (Real)strtod( value, NULL);
	if( fabs( *noise ) > 1 )
	{
		fprintf( stderr, "** invalid noise level=%g\n", *noise );
		return -1;
	}
	return 0;
}

/*******************************************************************************
 *
 *                  Program Helpers
 *
 ******************************************************************************/
#define MD_ERROR exit_value = -1; goto MD_END;


/*******************************************************************************
 *
 *                  Main Code
 *
 ******************************************************************************/
int main( int argc, char *argv[] )
{
	/***************************************************************************
	 *
	 * Variables Declaration
	 *
	 **************************************************************************/
	int             	exit_value     = 0;
	mds_frame       	frame          = MDS_FRAME_INIT;
	size_t          	num_particles  = 0;
	mds_step_type   	step_type      = MDS_STEP_VERLET;
	mds_step        	step           = MDS_STEP_INIT;
	mds_potential_type	potential_type = MDS_POTENTIAL_LJ;
	mds_mass_type       mass_type	   = MDS_MASS_ALEA;
	mds_units_type		units_type	   = MDS_UNITS_REAL;
	mds_potential_proc  potential_proc = mds_lj;
	void               *potential_args = NULL;
	int             	load_particles = 0;
	mds_verlet_data 	verlet_data;
	mds_beeman_data 	beeman_data;
	mds_dynamics    	dynamics        = { mds_dynamics_full, NULL };
	mds_dynamics        dynamics_simple = { mds_dynamics_full_simple, NULL };
	Real            	temperature    = 0;
	mds_vector_set  	position       = NULL;
	mds_vector_set  	velocity       = NULL;
	mds_lj_data     	lj;
	mds_lj_cutoff_data	lj_cutoff;
	mds_uvector     	cr;
	Real            	dt = 1.0e-3;
	Real            	t_run, t_old, t_max=5;
	Real            	t_therm = 0, K_therm = 0;
	register        	size_t i,iter;
	Real      			epsilon  = 0.995708;     
	Real      			sigma    = 0.3405;       
	Real                r_c      = MDS_F(2.5) * sigma;  /*-- current cutoff */
	mds_options     	options  = MDS_OPTIONS_INIT;
	mds_wtime           chrono;
	double              ellapsed = 0;
	Real                run_T, run_PV;
	Real                noise = 0;
	/*int 				tab_dens[];*/

	/***************************************************************************
	 *
	 * global initialization
	 *
	 **************************************************************************/
	mds_run_initialize();
	mds_wtime_init( &chrono );
	frame.pbc_for_r = mds_pbc_xyz;
	frame.pbc_for_f = mds_pbc_xyz;
	
	
	/***************************************************************************
	 *
	 * parse options
	 *
	 **************************************************************************/
	MDS_OPTIONS_DECL()
	{
		MDS_OPTION_NORMAL(f,         options,            "load options from file"),
		MDS_OPTION_FORCED(box,	     frame.box,   		 "simulation box"),
		MDS_OPTION_FORCED(np,   	 num_particles, 	 "number of particles"),
		MDS_OPTION_FORCED(dt,  		 dt,           		 "simulation time step"),
		MDS_OPTION_FORCED(tmax, 	 t_max,       	     "simulation max run time"),
		MDS_OPTION_NORMAL(mass, 	 mass_type,     	 "1|alea (default=massal)" ),
		MDS_OPTION_NORMAL(units, 	 units_type,     	 "real|adim (default=real)" ),		
		MDS_OPTION_NORMAL(step, 	 step_type,     	 "verlet|beeman (default=verlet)" ),
		MDS_OPTION_NORMAL(potential, potential_type,     "lj|lj-cutoff (default=lj)" ),
		MDS_OPTION_NORMAL(T,		 temperature,  	     "temperature (default 0)" ),
		MDS_OPTION_NORMAL(therm,	 t_therm,      		 "thermalization time (default 0 )" ),
		MDS_OPTION_NORMAL(noise,     noise,              "crystal noise level (default 0 )" )
	}
	MDS_OPTIONS_DONE();
	
	if( argc <= 1 )
	{
		MDS_OPTIONS_HELP(stderr);
		exit_value = 1;
		goto MD_END;
	}
	
	if( mds_options_parse( &options, argc, argv) < 0 )
	{
		mds_perror( stderr, "mds_options_parse()");
		MD_ERROR;
	}
	
	MDS_OPTIONS_DISPATCH(options,MD_ERROR);
	
	/***************************************************************************
	 *
	 * Check them...
	 *
	 **************************************************************************/
	
	
	fprintf( stderr, "md> num_particles= %u\n", (unsigned) num_particles );
	fprintf( stderr, "md> box=[ %g %g %g ]\n", frame.box.x, frame.box.y, frame.box.z );
	fprintf( stderr, "md> dt   = %g\n", dt );
	fprintf( stderr, "md> tmax = %g\n", t_max );
	

	/***************************************************************************
	 *
	 * prepare step functions
	 *
	 **************************************************************************/
	switch( step_type )
	{
		case MDS_STEP_VERLET:
			fprintf( stderr, "md> verlet step\n");
			mds_verlet_step_load( &step, &verlet_data );
			break;
			
		case MDS_STEP_BEEMAN:
			fprintf( stderr, "md> beeman step\n");
			mds_beeman_step_load( &step, &beeman_data );
			break;
	}
	assert( step.num_sets > 0 );

	/***************************************************************************
	 * 
	 * Choose units type
	 * 
	 **************************************************************************/
	
	switch( units_type )
	{
		case MDS_UNITS_REAL:
			fprintf( stderr, "md> units real\n");
			break;
			
		case MDS_UNITS_ADIM:
			fprintf( stderr, "md> units adim (warning m=1)\n");
			/* temperature = (MDS_R * temperature) / epsilon; */
			
			
			
			epsilon  = 1;
			sigma    = 1;
			MDS_R    = 1;
			mass_type = MDS_MASS_1;		/*warning automate choice for m if addim*/
			break;
	}	
	
	/***************************************************************************
	 *
	 * prepare potential
	 *
	 **************************************************************************/
	switch( potential_type )
	{
		case MDS_POTENTIAL_LJ:
			fprintf( stderr, "md> potential = lj\n");
			potential_proc = mds_lj;
			potential_args = &lj;
			mds_lj_init( &lj, epsilon, sigma );
			break;
			
		case MDS_POTENTIAL_LJCUTOFF:
			fprintf( stderr, "md> potential = lj-cutoff\n");
			potential_proc = mds_lj_cutoff;
			potential_args = &lj_cutoff;
			mds_lj_cutoff_init(&lj_cutoff, epsilon, sigma, r_c );
			
			break;
			
	}
	assert( step.num_sets > 0 );
		
	/***************************************************************************
	 *
	 * prepare frame and particles
	 *
	 **************************************************************************/
	
	if( load_particles )
	{
		/** TODO: load masses, positions, velocities... */
	}
	else 
	{
		/***********************************************************************
		 * adjust num_particles to crystal
		 **********************************************************************/
		num_particles = mds_crystal_close_packed_count( &cr, num_particles, & frame.box );
		fprintf( stderr, "md> adjusting to fcc: num_particles=%u\n", (unsigned) num_particles );
		if( mds_frame_load( &frame, num_particles, step.num_sets) < 0 )
		{
			mds_perror( stderr, "mds_frame_load( num_particles=%u, num_sets=%u )", (unsigned) num_particles, (unsigned) step.num_sets );
			MD_ERROR;
		}
		
	}
	assert( frame.num_sets >= 3 );
	assert( frame.sets != NULL  );
	
	position = frame.sets[ MDS_POSITION ];
	velocity = frame.sets[ MDS_VELOCITY ];
	
	if( ! load_particles )
	{
		/***********************************************************************
		 * build crystal
		 **********************************************************************/
		mds_crystal_build_fcc( position, frame.num_particles, &cr, & frame.box, noise);
		
		/***********************************************************************
		 * assign masses, here there are random.
		 **********************************************************************/
		switch( mass_type )
		{
			case MDS_MASS_1:
				fprintf( stderr, "md> mass = 1\n");
				for( i = 0; i < num_particles; ++i )
				{
					const Real mass           = 1;
					velocity[i].info.mass     = mass;
					position[i].info.inv_mass = MDS_F(1.0) / mass;
				}
				break;
				
			case MDS_MASS_ALEA:
				fprintf( stderr, "md> mass alea\n");
				for( i = 0; i < num_particles; ++i )
				{
					const Real mass           = 10 + 10 * mds_alea();
					velocity[i].info.mass     = mass;
					position[i].info.inv_mass = MDS_F(1.0) / mass;
				}
				break;
		}		

		/***********************************************************************
		 * assign gaussian velocities and frame.K
		 **********************************************************************/
		mds_frame_initialize_velocities( &frame, temperature);
	}
	
	
	/***************************************************************************
	 *
	 * Initialize time/accelerations
	 *void mds_density_calcul( mds_frame *frame );
	 **************************************************************************/
	iter  = 0;
	t_run = 0;
	t_old = t_run;
	K_therm = num_particles * ( MDS_F(1.5) * MDS_R * temperature);
	fprintf( stderr, "md> K_therm = %g\n", K_therm );
	mds_dynamics_call( &dynamics, &frame, potential_proc, potential_args);
	mds_step_init_dt( &step, &frame, dt);
	
	(void) mds_xyz_save( "md.xyz", position, frame.num_particles, "first", 0);
	(void) mds_file_save( "md.mds.gz", &frame, 0.0, 0 );
	
	{
		FILE *fp = fopen( "energy.txt", "wt" );
		if( fp )
		{
			mds_frame_NVE_T_and_PV( &frame, &run_T, &run_PV);
			fprintf( fp, "%g %.6e %.6e %.6e %.6e %.6e\n", t_run, frame.V, frame.K, frame.V + frame.K, run_T, run_PV);
			fclose( fp );
		}
	}
	MDS_MEMORY();
	mds_wtime_start( &chrono );
	while( t_run < t_max )
	{

		t_run = ++iter * dt;

		/*mds_step_predict( &step,  &frame ); 
		mds_dynamics_call( &dynamics, &frame, potential_proc, potential_args); */ /* V and W are updated */
		/*mds_step_correct( &step,  &frame );               */    /* K is updated here   */
		
		if( t_run <= t_therm )
		{
			frame.K = mds_velocity_rescale( velocity, num_particles, (Real)sqrt(K_therm/frame.K) );
		}

		if( (iter%100) != 0 )
		{
			mds_step_predict( &step,  &frame ); 
			mds_dynamics_call( &dynamics_simple, &frame, potential_proc, potential_args); 
			mds_step_correct( &step,  &frame ); 
		}
		/*if( t_run - t_old > 0.05 )*/
		else /*attention l'indice a 1 d'avance par rapport a t_run*/
		{
			mds_step_predict( &step,  &frame ); 
			mds_dynamics_call( &dynamics, &frame, potential_proc, potential_args);  /* V and W are updated */
			mds_step_correct( &step,  &frame );                                     /* K is updated here   */
			
			
			ellapsed += mds_wtime_query( &chrono );
			fprintf( stderr, "t= %8.4g |V=%.6e |K=%.6e |E=%.6e |<%8.4f> fps    \r", t_run, frame.V, frame.K, frame.K + frame.V , ((double)iter)/ellapsed );
			fflush( stderr );
			(void) mds_xyz_save( "md.xyz", position, frame.num_particles, "other", 1);
			(void) mds_file_save( "md.mds.gz", &frame, t_run, 1 );

			{
				FILE *fp = fopen( "energy.txt", "at" );
				if( fp )
				{
					mds_frame_NVE_T_and_PV( &frame, &run_T, &run_PV);
					fprintf( fp, "%g %.6e %.6e %.6e %.6e %.6e\n", t_run, frame.V, frame.K, frame.V + frame.K, run_T, run_PV);
					fclose( fp );
				}
			}
			/*t_old = t_run;*/
			mds_wtime_start( &chrono );
		}
	}
	fprintf(stderr, "\n\nmd> done.\n");
	
	mds_nombre_accumulation( &frame, sigma );
	/*mds_nombre_accumulation_radial( &frame );*/

	/***************************************************************************
	 *
	 * End of the program
	 *
	 **************************************************************************/
	goto MD_END;
MD_END:
	
	/** release resources */
	mds_options_quit( &options );
	mds_frame_quit( &frame );
	
	MDS_MEMORY();
	return exit_value;
}
